#Get objects
old_ws = ls()

#############
#Tables NVMS#
#############

useTable_nvms = mergedTable_nvms[!is.na(NVM_Violence_1_count) & (sample_ni_it)]

useTable_nvms_2004 = useTable_nvms[election_cycle %in% 2004]
useTable_nvms_2009 = useTable_nvms[election_cycle %in% 2009]
#################
#Events Analysis#
#################

#Model specifications
nb_specs = expand.grid(x = 'IT_win', y = c('NVM_Violence_1_count', 'NVM_Violence_1_deaths'), stringsAsFactors = F)
lm_specs = expand.grid(x = 'IT_win', y = c('NVM_Violence_1_binary', 'NVM_Violence_1_deaths_binary'), stringsAsFactors = F)

#Estimate models
nvms_nb_models_2004 = Map(function(y,x) glmnb_loop(y,x, useTable_nvms_2004),
                     nb_specs$y, nb_specs$x)

nvms_lm_models_2004 = Map(function(y,x) lm_loop(y,x,  useTable_nvms_2004),
                     lm_specs$y, lm_specs$x)

nvms_nb_models_2009 = Map(function(y,x) glmnb_loop(y,x, useTable_nvms_2009),
                          nb_specs$y, nb_specs$x)

nvms_lm_models_2009 = Map(function(y,x) lm_loop(y,x,  useTable_nvms_2009),
                          lm_specs$y, lm_specs$x)


nvms_list_2004 = c(nvms_nb_models_2004, nvms_lm_models_2004)
nvms_list_2009 = c(nvms_nb_models_2009, nvms_lm_models_2009)

nvms_se_2004 = lapply(nvms_list_2004, function(x) cluster_errors(x, useTable_nvms_2004[, cluster]) %>% diag %>% sqrt)
nvms_se_2009 = lapply(nvms_list_2009, function(x) cluster_errors(x, useTable_nvms_2009[, cluster]) %>% diag %>% sqrt)

########################
#Tables Violence: PODES#
########################

useTable_podes = mergedTable_podes[!is.na(Violence_1_count) & (sample_ni_it)]

useTable_podes_2004 = useTable_podes[election_cycle %in% 2004]
useTable_podes_2009 = useTable_podes[election_cycle %in% 2009]
#############
#Full sample#
#############

#Model specifications
nb_specs = expand.grid(x = 'IT_win', 
                       y = c(paste0('Violence_',1,'_count'),
                             paste0('Violence_',1,'_deaths')), stringsAsFactors = F
)

lm_specs = expand.grid(x = 'IT_win', 
                       y = c(paste0('Violence_',1,'_binary'), 
                             paste0('Violence_',1,'_deaths_binary')), stringsAsFactors = F
)

#Estimate models
podes_nb_models_2004 = Map(function(y,x) glmnb_loop(y,x, useTable_podes_2004),
                      nb_specs$y, nb_specs$x)

podes_lm_models_2004 = Map(function(y,x) lm_loop(y,x,  useTable_podes_2004),
                      lm_specs$y, lm_specs$x)


podes_nb_models_2009 = Map(function(y,x) glmnb_loop(y,x, useTable_podes_2009),
                      nb_specs$y, nb_specs$x)

podes_lm_models_2009 = Map(function(y,x) lm_loop(y,x,  useTable_podes_2009),
                      lm_specs$y, lm_specs$x)


podes_list_2004 = c(podes_nb_models_2004, podes_lm_models_2004)
podes_list_2009 = c(podes_nb_models_2009, podes_lm_models_2009)

podes_se_2004 = lapply(podes_list_2004, function(x) cluster_errors(x, useTable_podes_2004[, cluster]) %>% diag %>% sqrt)
podes_se_2009 = lapply(podes_list_2009, function(x) cluster_errors(x, useTable_podes_2009[, cluster]) %>% diag %>% sqrt)


############
#Make Table#
############

table_list_2004 = c(podes_list_2004, nvms_list_2004 )
table_se_2004 = c(podes_se_2004, nvms_se_2004)

table_list_2009 = c(podes_list_2009, nvms_list_2009)
table_se_2009 = c(podes_se_2009, nvms_se_2009)

#create table 2004
n_clusters = useTable_podes_2004[(nvms_sample), cluster] %>% unique %>% length
n_clusters_2 = useTable_podes_2004[, cluster] %>% unique %>% length

note_text = paste("Count models use negative binomial regression; binary outcomes use OLS. Standard errors are clustered by constituency-clusters. In the NVMS and PODES samples, there are ", n_clusters, "and ", n_clusters_2, " clusters, respectively. Observations are constituencies in which the last seat was contested by Islamist and secular nationalist parties with a margin less than 1 percent.")

table = stargazer(table_list_2004, type = 'latex', 
                  title = "Estimated Effects of Islamist Victory on Religious Violence (2004 Only)",
                  label = 'tab:main_violence_ni_it_2004',
                  model.names = F,
                  model.numbers = T,
                  column.separate = c(4,4),
                  column.labels = c('PODES', 'NVMS'),
                  multicolumn = T,
                  dep.var.labels = rep(c("Events","Deaths"), 4), 
                  add.lines = list(c('Count', rep(c('Y', 'Y', 'N', 'N'), 2)),
                                   c('Binary', rep(c('N', 'N', 'Y', 'Y'), 2))),
                  covariate.labels = c("Islamist Win"),
                  star.cutoffs = c(0.1, 0.05, 0.01),
                  #float.env = 'sidewaystable',
                  notes.align = 'l',
                  font.size = 'scriptsize',
                  keep.stat = 'n',
                  style = 'apsr')

write_latex(table[-10] %>% append(table[10], after = 10) %>% append("\\cmidrule(lr){2-5}\\cmidrule(lr){6-9}", after = 10), note_text, './output/tables/table_e8.tex')

#create table 2009
n_clusters = useTable_podes_2009[(nvms_sample), cluster] %>% unique %>% length
n_clusters_2 = useTable_podes_2009[, cluster] %>% unique %>% length

note_text = paste("Count models use negative binomial regression; binary outcomes use OLS. Standard errors are clustered by constituency-clusters. In the NVMS and PODES samples, there are ", n_clusters, "and ", n_clusters_2, " clusters, respectively. Observations are constituencies in which the last seat was contested by Islamist and secular nationalist parties with a margin less than 1 percent.")

table = stargazer(table_list_2009, type = 'latex', 
                  title = "Estimated Effects of Islamist Victory on Religious Violence (2009 Only)",
                  label = 'tab:main_violence_ni_it_2009',
                  model.names = F,
                  model.numbers = T,
                  column.separate = c(4,4),
                  column.labels = c('PODES', 'NVMS'),
                  multicolumn = T,
                  dep.var.labels = rep(c("Events","Deaths"), 4), 
                  add.lines = list(c('Count', rep(c('Y', 'Y', 'N', 'N'), 2)),
                                   c('Binary', rep(c('N', 'N', 'Y', 'Y'), 2))),
                  covariate.labels = c("Islamist Win"),
                  star.cutoffs = c(0.1, 0.05, 0.01),
                  #float.env = 'sidewaystable',
                  notes.align = 'l',
                  font.size = 'scriptsize',
                  keep.stat = 'n',
                  style = 'apsr')

write_latex(table[-10] %>% append(table[10], after = 10) %>% append("\\cmidrule(lr){2-5}\\cmidrule(lr){6-9}", after = 10), note_text, './output/tables/table_e9.tex')



#Clean up
drop = setdiff(ls(), c(old_ws)) 
rm(list = drop)
